CU-TP-1198, INT-PUB-1 1-024, RIKEN-QHP-3 

Lattice Monte Carlo calculations for unitary fermions in a 

harmonic trap 

Michael G. Endreffl 

Physics Department, Columbia University, New York, NY 10027, USA and 
Theoretical Research Division, RIKEN Nishina Center, Wako, Saitama 351-0198, Japan 

David B. Kaplanj^ Jong- Wan Leej^ and Amy N. Nicholsonl^ 

Institute for Nuclear Theory, University of Washington, Seattle, WA 98195-1550, USA 

(Dated: November 4, 2011) 

Abstract 

We present a new lattice Monte Carlo approach developed for studying large numbers of strongly 
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harmonic trap. In place of importance sampling, our approach makes use of high statistics, an 
improved action, and recently proposed statistical techniques. We show how improvement of the 
lattice action can remove discretization and finite volume errors systematically. For = 3 unitary 
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benchmarks we reproduce precision calculations of = 3,..., 6 unitary fermions in a harmonic 
trap to within our ~ 1% uncertainty. We then use this action to determine the ground state 
energies of up to 70 unpolarized fermions trapped in a harmonic potential on a lattice as large as 
64^ X 72. In contrast to variational calculations we find evidence for persistent deviations from the 
thermodynamic limit for the range of N considered. 
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I. INTRODUCTION 



Developing a predictive understanding of strongly interacting many-body systems is one 
of the most difficult and potentially rewarding challenges in physics. A paradigm for this 
problem in perhaps its purest form is to determine the behavior of a gas of unitary fermions 
(for a brief overview, see [T]). These are nonrelativistic fermions with zero range interactions 
tuned such that the two-body s-wave scattering length diverges. Thus the s-wave phase shift 
satisfies 6{k) = tt/2 for all k and the field theory describing the many-body system is at 
a conformal fixed point^; in 1998 it was suggested that unitary fermions could serve as 
the starting point for an effective field theory expansion for nuclear physics |21 E]- Since 
then the unitary fermion gas has been created and studied experimentally by trapping 
atoms tuned to a Feshbach resonance by means of an applied magnetic field, exhibiting 
collective effects interpolating between the well understood phenomena of BCS pairing and 
Bose-Einstein condensation [1HT3]. The nonperturbative nature of the strongly coupled 
interaction between unitary fermions poses a nontrivial challenge for theory, and numerical 
simulation has played an essential role in making progress. A large body of recent theoretical 
work exists for unitary fermions, both analytical [T1H20] and numerical |211H3] . 

In this paper we describe a new lattice approach for simulating unitary fermions, and 
determine the ground state energies for up to 70 unitary fermions in a harmonic trap on 
lattices as large as 64^ x 72, allowing for an extrapolation to the infinite volume limit. 
This significantly extends preliminary findings published in the Lattice 2010 conference 
proceedings |H^fl6] . building on the lattice construction of |17]. In addition, for this work 
we have made several improvements, including the use of a Galilean-invariant interaction ^ 
for tuning to unitarity and reducing time discretization errors in the implementation of the 
harmonic oscillator; these are outlined in Sec. |TT} 

Our approach differs from previous numerical studies in several ways: 

^ Since the underlying theory is conformal, at nonzero chemical potential /i and h — 1, all dimensionful 
quantities, such as the ground-state energy and pairing gap A, are given as pure numbers times the 
function of fi and the fermion mass M combined to give the corresponding dimension. 

^ By Galilean invariant, we mean that the interaction is only a function of the transferred three-momentum 
between interacting particles, although that momentum is necessarily discrete and periodic on the fi- 
nite volume lattice. Momentum dependent separable interactions, for example, would not be Galilean 
invariant. 
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• The theory is defined on a four dimensional Euchdian lattice, and fermion-fermion 
interactions are induced by an auxiliary scalar field (f). We compute A^-fermion corre- 
lators in the background field, then average observables over an ensemble of these 
fields - in much the same way one computes the hadron spectrum in lattice QCD. 
Unlike some approaches [HHlS], our computation is not variational in nature, and so 
our result for the ground state energy does not depend on an accurate parametrization 
of the many-body ground state wavefunction. In practice, however, using good sources 
and sinks for the correlators is necessary to achieve this goal, blurring the boundary 
between unconstrained and variational calculations when N is large. 

• We formulate the lattice action in such a way that the fermion determinant is indepen- 
dent of the auxiliary field so that the so-called "quenched approximation" is exact, 
greatly simplifying the computation. This requires open boundary conditions in the 
temporal direction, and we can therefore only study properties at zero temperature. 

• We do not use importance sampling (that is, we do not include the correlator we are 
trying to compute as part of the measure for 0). Instead, our ensemble consists 
of random Z2 valued variables living on the time-like lattice links, and therefore is 
extremely cheap to generate (see [H] for a detailed discussion of the scaling of our 
algorithm with volume and number of fermions). The price we pay is that we face 
a serious distribution overlap problem that cannot be overcome simply by increasing 
statistics ^. 

• The distribution overlap problem is identified as arising from heavy-tailed distributions 
for our correlators, similar to what is seen for conductance electrons in disordered me- 
dia near the Anderson localization transition. We have developed a statistical method 
for greatly ameliorating the problem, as discussed in a separate paper, Ref. [i8] . 

• We use a greatly improved lattice action that exactly reproduces single particle disper- 
sion relations up to a momentum cutoff related to the inverse lattice spacing as well as 

Due to an unfortunate choice of nomenclature, the "overlap problem" commonly refers to one of two 
unrelated problems, both of which concern us here. The first is the poor overlap between the true ground 
state and the choice of interpolating operators, whereas the latter is the poor overlap between the path- 
integral probability measure and the dominant part of the operator being estimated. We will refer to the 
former as a "interpolating operator overlap problem" and the second as a "distribution overlap problem." 
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the first several two-particle energy levels in a box with zero lattice spacing. We show 
that the volume dependence we find for the energies of two-body states are consistent 
with fermions having the first four or five terms in the effective range expansion tuned 
to zero. Thus our fermions are much closer to the unitarity limit than have ever been 
studied before for > 3 particles, and as a result we have small discretization errors 
and do not have to extrapolate our results to zero range, as do most simulations. 

We have formulated this theory both for unitary fermions in a box ( "untrapped" ) or in a 
harmonic potential ("trapped"). In this paper we will present only the results for trapped 
fermions, leaving the untrapped results for future publication |19], although we use results 
for two and three untrapped fermions to help establish the validity of our method. 

The organization of this paper is as follows. In Sec. [IT] we describe the theoretical details 
of our lattice construction, including notational conventions, lattice parameter tuning meth- 



ods, and an analysis of discretization errors. In Sec. Ill we present ensemble details and 



measurement results for the ground state energies of up to 70 unpolarized unitary fermions 



confined to a harmonic trap. We conclude in Sec. IV with a summary of results and a 
discussion of possible future applications of our lattice construction. More technical details 
are provided in appendices: Appendix |X] gives details about tuning the lattice interaction; 
Appendix [B] describes how we construct our multi-fermion correlators which incorporate 
pairing correlations; Appendix [O explains our strategy for extracting accurate estimates of 



the multi-fermion energies using cumulant expansion techniques of Ref. |18]; Appendix D 
provides details of our simulation, including various numerical checks performed in order to 
verify the correctness of our code. 



II. LATTICE CONSTRUCTION 



A. Action, notation and conventions 



The starting point for our construction is a highly improved variant of the nonrelativistic 
Euclidean-time lattice action proposed in [47j: 



S = Khl J2 



1 , 



(1) 
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This action describes two species of one- component interacting fermions = (V'^, ip'^) with 
equal mass M defined on a T x L'^ lattice, with the temporal and spatial lattice spacings 
given by bj. and bs, respectively. For convenience, we work primarily in lattice units, where 
bg = br = 1, however in some sections we restore the lattice spacings in order to discuss 
temporal and spatial discretization errors. Throughout this work, we consider a lattice with 
open boundary conditions in the time direction with time labeled by integers r G [0, T — 1], 
and periodic boundary conditions in the spatial directions with position labeled by integers 
Xj G [— L/2, L/2 — 1], for j = 1,2, 3. As a result of using open temporal boundary conditions, 
the utility of our lattice action is limited to studies at zero temperature. In addition, this 
choice of boundary conditions forbids the introduction of a chemical potential and we work 
in the canonical, rather than grand-canonical ensemble. 

The derivative operator dr appearing in Eq. [ijrepresents a backward difference operator in 
time, i.e., (9t-'?/')x,t = "^x.t — "^x.t-i, whereas represents a lattice gradient operator defined 
so as to give a perfect continuum-like single particle dispersion relation for free fermions. 
This kinetic term is highly nonlocal, although as will be described below, the nonlocality 
poses no challenge in a numerical simulation of Eq. [T} 

A four-fermion contact interaction is achieved via the introduction of a stochastic auxil- 
iary scalar field 0x,t associated with the time-like links of the lattice. This field is chosen to 
satisfy the conditions 



where the expectation value represents ensemble averaging over 0, and in this work the (j) 
distribution is taken to either be unit- variance Gaussian or The point-split character of 
the interaction ensures that scattering propagates fermions forward in time by one unit. This 
choice, along with the absence of fermion propagation in the negative temporal direction 
and open boundary conditions in time, ensures that no closed fermion loop depends on 0. 
A consequence is that the fermion determinant is ^-independent and has no effect on the 
measure for 0, greatly simplifying numerical simulation of Eq. [1} 

The operator Cxx' = C(x — x') acts only in space and is taken to be real, symmetric, 
local, and invariant under lattice translations; it can be thought of as a differential operator 
acting on which allows the interaction between fermions induced by exchange to depend 
on the transfer momentum. Not only does this give us a momentum-dependent interaction 





(2) 
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we can tune to attain unitarity, but it is also Galilean invariant in that it depends only on 
the difference between the ingoing and outgoing fermion momenta. This is important, since 
tuning a non-Galilean invariant interaction to give unitarity in one frame would lead to non- 
unitary fermions in another, and boosted pairs of particles would see an interaction which 
did not correspond to unitarity. Integrating out the auxiliary field yields the four-fermion 
interaction 

where ('?/'V')x,r = 'V^x,rV'x,r-i, and we have used the Hermiticity of C . We may express Eq. [l] 
succinctly as S* = tpKip, where the time components of the fermion matrix K are given in 
block-matrix form by: 



K = 



with 
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(4) 



Note that the x matrices D, X, C and $(r) act only in space and that $(r) is a 
diagonal matrix with statistically independent random elements 0x(t). 

We choose to realize the lattice Laplacian in such a way that D has the following form 
in momentum space (|27]): 

f eP'/(2^'^) IpI < A 
^PP' = f^p.p' X < , (6) 

I oo IpI > A 

where pj = 2Timj/L for integers mj G [— L/2,L/2 — 1] and j = 1,2,3. The parameter 
A = 7r X (1 — 10~^) is a hard momentum cutoff imposed on the fermions; a small shift away 
from TT has been introduced in the cutoff in order to avoid inclusion of momenta lying on 
the very edge of the Brillouin zone (BZ). For free fermions, X = 1 and the propagator is 
just a transfer matrix, which in momentum space has the form 

Ke(o,-)]pp, = [^-^pp, - We-^(^^^^(A - IpI) (7) 
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and yields the exact one-particle energy, E{p) = p^/2M. So we see that the choice Eq. 
(|6| is designed to give the exact one-particle dispersion relation up to a momentum cutoff 
IpI = A, beyond which the fermions do not propagate. Imposing the A cutoff just within the 
Brillouin zone boundary was necessary to reconcile the exact continuum dispersion relation 
with the periodicity of the reciprocal lattice. 

For the interaction we take in momentum space 

Cpp' = 5p,p' X ^ , (8) 

[ C(A) IpI > A 

where below A, C(p) is an analytic function of p^ which we adjust to construct the desired 
continuum phase shift for two-particle scattering (for example, the constant 5 = 7r/2 phase 



shift for unitary fermions). How we tune C is discussed in Sec. II B and Sec. II C 

In order to simulate the partition function defined by Eq. [1} it is necessary to first integrate 
out the fermionic degrees of freedom, yielding an effective action involving only the auxiliary 
field. The resulting partition function is given by 



j [c/0]p(0) det K , (9) 



where 



e 2<Px Gaussian 

The corresponding expectation value of an arbitrary operator 0{ip,il!) is given by: 

m, n = ^ j [#]P(0) det K d{K-') , (11) 

where OlK"^) is some new calculable operator that depends implicitly on (p through the 
propagator K^^. Both O and O may have explicit dependence on as well. Since K is an 
upper triangular block matrix, its determinant is given by the product of determinants of its 
diagonal blocks, det K = (det -D)"^, which is independent of the auxiliary field. Therefore the 
full numerical simulation of the partition function with action given in Eq. [T] is equivalent 
to a quenched simulation, with expectation values given by: 

{0{^,ij)) = / mp{<j>)d{K-') , (12) 

■^quenched J 
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where ^quenched = J[d4>]p{'P) IS the quenched partition function. Note that the absence of a 
nontrivial probabihty measure for the auxihary field ensures that the path integral is free of 
the sign problem. 

Because K is upper triangular in form, interacting fermion propagators measured from 
time slice zero to time slice r may be expressed exactly as a sequence of applications of 
and X operators, resulting in a simple recursive formula: 

K-\t; 0) = D'^X{t - 1)K-\t - 1; 0) , (13) 

with i^~^(0;0) = D~^. The form of this result is evident from the fact that there are no 
time-like closed fermion loops, which is a consequence of using open boundary conditions and 
from the absence of anti-particles in the nonrelativistic theory. Inversion of the nonlocal D 
operator and application of the X{t) operator may be performed efficiently with fast Fourier 
transforms (FFTs); it is this feature that allows us to use the perfect dispersion relation and 
momentum dependent interaction defined in Eq. |6] and Eq. |8j 

B. Transfer matrix formalism 

Multi-fermion correlation functions C(r) are obtained from an ensemble average of direct 
products of propagators 

)C-\t- 0) = K~\t- 0) ® . . . ® K-\t] 0) , (14) 

^ V ' 

N 

which are sandwiched between properly antisymmetrized A^-fermion initial and final states 
(i.e., interpolating fields associated with time slices zero and r, respectively). We will refer 
to the initial and final states as sources and sinks, respectively. We may translate our 
lattice action in Eq. [T] into Hamiltonian language by noting that the expectation value of 
/C^^(r; 0) is just the Euclidean time evolution operator for a system of particles. Since the 
single particle propagator is itself a product of uncorrelated random matrices (because 
the auxiliary field probability measure is separable in time), the multi-fermion correlation 
function will factor into a matrix product of ensemble averages. If we define the matrix: 

T = V-^/\l-V)V~^^\ (15) 
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where 

V = , (16) 

N 

and 

(1-V) = (X(T)(g)...(g)X(r)) , for every r (17) 

^ V ' 

N 

are dimensional matrices, then the A^-fermion correlator may be written in the highly 
suggestive form: 

(/C-i(r; 0)) = V^^^^ryV-^^^ , (18) 

and we may identify T as a transfer matrix and H = — In T as a Hamiltonian for the 
iV-fermion system, provided T is Hermitian and positive^. 

A general expression for the multi-particle interaction V may be computed analytically 
from Eq. [5] and Eq. [T7| by explicit integration of the auxiliary fields. The expression is 
somewhat complicated for large numbers of particles and will therefore not be explicitly 
derived here. Observe, however, that although the auxiliary field interaction X{t) involves 
a square root of the operator C, the multi-particle interaction V is in fact an analytic 
function of momenta. This is due to the presence of momentum conserving delta functions 
which ensure that y/C always comes in pairs; in terms of Feynman diagrams, there are 
identical factors of y/C at each end of the propagator, only depending on the magnitude 
of the momentum flowing through that propagator. This property is generally true for 
any A^-particle system since only an even number of insertions of the interaction survive 
integration over the auxiliary fields; it is also evident from the right-hand-side of Eq. [3] 



In the case of two fermions, where Ni = = 1, the transfer matrix defined by Eq. 15 
may be evaluated in momentum space and is given by: 

\q q I / IP p ; - g(q^2^qt2+piVpt^)/{4M) ' ^^^^ 

for momenta below the cutoff A. C(p) is a periodic function of the operator p for |p| < A 
which we choose to expand in a convenient basis of local functions: 

'^(P) = C2n02n{p) , (20) 

n=0 



^ This is a stronger condition than necessary; if T is hermitian but not positive then TT = T^T 'is Hermitian 
and positive, guaranteeing that a sensible definition of the Hamiltonian will exist with a time step of 267- . 
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with unknown coefficients C2n to be determined from scattering data. Our choice of basis 
functions is: 

f f 1 - e-p'/^^o V |p|<A 
0,^{p)=M-x{) J (21) 

[(l-e-^'/^«) |p|>A 

for p within the first Brillouin zone, and periodic from one Brillouin zone to the next. The 
basis functions behave as 02n{p) ~ P^" for small p^ <ti Mq and tend to a constant for 
p^ > Mo; this basis was chosen to approximate continuum 2-body contact interactions with 
2n derivatives for low transfer momentum, while not getting excessively big for momenta 
at the edge of the Brillouin zone. Throughout this work we take Mq = M, and both to be 
0(1) in lattice units. 

In the special case where Nq = 1 the only operator in the sum Eq. |20] is Oq which 
is constant, and the two-fermion transfer matrix may be diagonalized analytically on the 



finite volume lattice. All nonzero total momentum eigenstates of Eq. 19 correspond to plane 
waves, whereas the zero total momentum eigenstates are given by 

(pVI^^) oc ^_^^^^./^_^ Wpt,o (22) 

where p = |p^| = |p^|. The corresponding energy eigenvalues E/^ are given by solutions to 
the integral equation 



= (23) 

p<A 



^0 



which, for every value of p"^, admits a single bound state for any value of Co > at finite 
volume. This negative energy state becomes a scattering state in the infinite volume limit 
for < Co < Ccrit and a bound state for Ccrit < Co, where Ccrit is an M-dependent critical 
value; tuning Co — )■ Ccrit yields a zero energy bound state at infinite volume, corresponding 
to unitarity and the continuum limit of the lattice theory. 

In the case where Nq > 1, even semi- analytic solutions for the C2„ coefficients are not 



feasible, but they may be determined numerically by explicit diagonalization of Eq. 19 It 
is helpful to restrict the transfer matrix to the zero center-of-momentum subspace, thus 
reducing the dimensionality of the matrix from down to a more manageable size of L^. 



A further reduction in the dimensionality of Eq. 19 may be achieved by projecting the zero 



center-of-momentum part of the transfer matrix onto appropriate representations of the 
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octahedral group Oh (e.g., in the case of s-wave scattering, the trivial representation Af). 
Performing such a projection makes numerical diagonalization feasible for lattices at least 
as large as L = 64, which is the maximum lattice size we consider in our numerical studies. 

C. Parameter tuning 

Unitary fermions in the continuum are a conformal system, while a lattice simulation 
necessarily involves finite lattice spacing and volume, both breaking conformal symmetry. 
Critical to a numerical simulation is the ability to tune the interactions to unitarity and 
control the systematic errors. In contrast to chiral symmetry in lattice QCD, for example, 
there is no phase transition associated with unitarity, despite the enhanced symmetry, and 
so there is no general feature in the A^-body spectrum that allows one to easily evaluate how 
far one is from unitarity. It is important therefore to collect as many results as possible about 
unitary fermions in the continuum that are known exactly or to high numerical precision in 
order to facilitate the tuning of the lattice action and to control systematic errors. 

What is known exactly about unitary fermions in the continuum is (i) the spectrum of 
two unitary fermions in a box of size L [501453] : (ii) the spectrum of two and three unitary 
fermions in a harmonic trap [18J; (iii) the scaling dimension of local composite operators 
involving unitary fermions^. Not known exactly but determined to high numerical accuracy 
are (iv) the few lowest energy levels for three unitary fermions in a box, extrapolated from 
a lattice Hamiltonian diagonalization very close to the continuum limit, with lattice size up 
to L = 50 |29]; and (v) the ground state energies for 4, 5, 6 unitary fermions in a harmonic 
trap, obtained by solving the Schrodinger equation [57]. The ground state energy for = 4 
fermions in a box has also recently been precisely studied by several methods in Ref. [37], 
but involves extrapolation to the continuum from very small lattices, L < 8, which makes 
the evaluation of potential systematic errors difficult. 

Our strategy for utilizing this information to tune our lattice action and estimate the size 
of systematic errors is to adjust our C2n coefficients to correctly reproduce the low-lying two- 
particle spectrum in a box in the continuum, subsequently showing that we can reproduce 

^ The scaling of two-body operators was determined in Ref. [2l|3] (see also [54]); the scaling of low dimension 
three-body operators was first analyzed by Griesshammer |55l 156) , and a beautiful general analysis was 
subsequently supplied by Nishida and Son {17) . 
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the correct volume scaling relations of measured energies, as well as the precisely known 
ground state energies for 3-fermions in a box or 3-6 trapped fermions. Here we discuss 
the tuning and energy levels of two and three untrapped fermions; our results for few-body 



trapped fermions are discussed in Sec. Ill 



1. Tuning and scaling of low-lying 2-body untrapped energy levels 

The two-particle energies E for s-wave particle pairs in a box with zero net momentum 
and phase shift 60 are given by the solutions to 



p cot 60 = —S{ri) , S{ri) = lim 

TtL A-s>oo 



|j|<A 



4:71 A 



(24) 



where j is an integer three- vector, rj = (pL/27r)^, and p is related to the energy hj E = p'^/M 
[SOHSS]- If scattering is due to short range interactions, then pcot^o is analytic in at 
sufficiently low p and one has the effective range expansion, 

pcot6o = -^ + ^rop^ + rip^ . . . , (25) 

where a is the scattering length, rg is the effective range, and ri, with dimension of volume. 



is what we will call the shape parameter. By means of Eq. 24, knowledge of the energy 
eigenvalues for the low-lying two-particle modes in a box can be used to determine effective 
range expansion parameters. Conversely, given a target set of effective range expansion 



parameters, we can tune our operator coefficients in Eq- (20) of our lattice theory 
until we attain the correct low-lying energy eigenvalues. This general tuning procedure was 
introduced in [5B] ■ For unitary fermions in the continuum we set p cot (5o = on the lefthand 
side of Eq. ( [24) ) and find the solutions rjl to the equation S{r]l) = 0. The function S{t]) is 
shown in Fig. [T| and the roots rjl correspond to the points where the function crosses the r) 
axis. The first 27 solutions are listed in Table [l]^. 

On the lattice the energy eigenvalues are defined from A = e^''^^, where A are the eigenval- 
ues of the two-particle transfer matrix discussed above and br is the temporal lattice spacing. 



^ To compute the ry^ it is very helpful to recognize that the number of integer three vectors j with equal 
norm is given by the coefhcient of x'^' in the Taylor expansion of [^3(0, a;)] where 63{u,x) is one of the 
Jacobi theta functions. 



12 



150 
100 
50 

S 

-50 
-100 
-150 



-4 -2 2 4 6 8 10 

FIG. 1: A plot of the three-dimensional (^-function Sijj). 
TABLE I: First 27 roots r/* (A; = 1, . . . , 27) of 5(7?). 
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-0.0959007 
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7.1962633 
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15.3537376 
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23.0194729 
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0.4728943 


9 


8.2879537 


16 


16.1218254 


23 


24.3306210 
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1.4415913 


10 


9.5345315 


17 


17.5325416 


24 


25.3016129 
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2.6270076 


11 


10.5505341 


18 


18.6053932 
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26.6803601 
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3.5366200 


12 


11.7014958 


19 


19.5186394 


26 


27.8780020 
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4.2517060 


13 


12.3102392 


20 


20.4033187 


27 


29.6156511 
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5.5377008 


14 


13.3831152 


21 


21.6944179 







Spatial discretization effects make it impossible to exactly reproduce the continuum rjl on 
the lattice. For one thing, there are an infinite number of rjl while the lattice transfer matrix 
has only a finite number of eigenvalues. Furthermore, since the lattice restricts how easily 
fermions can get close to each other — effectively creating a repulsive interaction — the 
phase shift for lattice unitary fermions necessarily falls below tt/2 for large lattice momenta, 



and pcot Sq as computed from Eq. (24) gets large. So the best one can do is tune a number 
No of the C2n coefficients to reproduce the lowest No solutions r]l. Details of how this 
tuning was performed numerically are provided in Appendix |X| In Table [TT| we give as an 
example the results for tuning operators for an L = 32 lattice with mass M = 5. 

Once we have tuned the operator coefficients, we can compute all eigenvalues of the 2- 



particle transfer matrix relevant for continuum s-wave scattering and use Eq. 24 to determine 
pcot 6o. Fig. [2] shows the result of this exercise for the successive tunings of Table [Tij In the 
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TABLE II: Results for tuning Nq C2n coefficients for an L = 32, M = 5 lattice. Uncertainties in 
the coefficients reflect a numerical uncertainty in 77^ at 0(10" 7). 



No 


Co 


C2 


C4 




1 


0.6815346(1) 








2 


0.466516(2) 


0.0856007(8) 






3 


0.489085(8) 


0.00853(2) 


0.020778(6) 




4 


0.50142(5) 


0.00958(3) 


0.00350(8) 


0.00430(2) 



left panel we show that p cot 6q over a wide range of momenta, extending well beyond 
that of the < 4 lowest eigenvalues we used to tune the C2n- 

Having p cot 5q look progressively flatter with each tuning is only a qualitative indication 
that we are attaining unitarity with improvement at each order. It is not advisable to try to 
fit this curve with a polynomial to extract effective range expansion coefficients; the reason 
is that the lattice function is only defined at discrete points, and one expects a finite - but 
unknown - radius of convergence for the effective range expansion. As a result it is possible 
to extract wildly different effective range coefficients from a polynomial fit, depending on 
the order of the fit and its momentum range. The situation is clarified in the right panel 
of Fig. |2] which plots pcot^o on a In- In plot. This plot shows clear evidence that with each 
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FIG. 2: Left: pcot5o as computed from exact lattice 2-particle energy eigenvalues using Liischer's 
formula, with the first Nq terms in the effective range expansion tuned to zero for Nq = 1, • • • ,4. 
Right: Same data on a In-ln plot along with expected r/ scaling (dashed lines) for various Nq. Data 
is from an L = 32 and M = 5 lattice. 
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FIG. 3: Succesful tuning of effective range parameters may be seen in the L dependence of indi- 
vidual energy eigenvalues for two particles in a box. Here we see agreement with Eq. [28] for the 
L-dependence (in lattice units) of levels rj^ and ryg which were not tuned. 

successive tuning we are setting successive terms in the effective range expansion to zero. 
Furthermore, the convergence of the dashed lines in the plot at r/ ~ 30 demonstrates that 
the radius of convergence for the effective range expansion is ?7 ~ 30, with deviations of the 
plotted points from the dashed lines indicating significant breakdown of the expansion at 
f] > 15, or \p\ ~ 0.76/bs. Note that for free fermions, rj is an integer that denotes the energy 
shell, and that a degenerate fermi gas filled to the ?7 = 15 shell would contain 251 fermions 
of each spin, far above the number of fermions we actually are able to study ^. 

Another way to see if the tuning procedure is successful is to look at the L-dependence 
of the low-lying energy eigenmodes on the lattice. Assume that we have tuned p cot 6q so 
that the leading term in the effective range expansion is 

where r„_i has dimensions (length)^""^, and that rjk are the solutions to S{ri) = nLpcot6o, 
while as before, the rjl are the unitary limit solutions to S{r]) = 0. For sufficiently small 
rjk — rjl we have S{rik) — Ck{rjk — rjl) where Ck are the slopes of S where it intersects the 
r^-axis in Fig. [T} Thus we find 

i(27r)2"+iLi-2V„_i(r^*)" ^ - vl) (27) 

^ The scattered behavior of the lowest 77 points in the right panel of Fig. [2] seem to indicate the difhculties 
with our procedure when we attempt to tune too many C2n parameters. 
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or 



(28) 

Thus the prediction is that a plot of L — 1^ should scale like L"*-^""^) when n terms in 
the effective range expansion have been tuned away. Note that because of the L~^" factor in 
the above equation, the effects of a small residual term at lower order in the effective range 
expansion will dominate at sufficiently large L. We have computed the low-lying energy 
eigenvalues for two particles on lattices of a number of different sizes, and in Fig. [3] we plot 
the results for energy levels 775 and 779, both at higher shells than were used in our tuning 



procedure. The scaling of Eq. [28] is evident in these plots: at each successive tuning we see 
that the L dependence is steepened by an additional factor of L^^. An interesting exception 
is for ?75 with four parameters tuned and L > 22; there we see points flattening out to 
perhaps an slope, suggesting that a small residual shape parameter ri is beginning to 
dominate at that point. We can use this deviation, Eq. (28), the value of rfl from Table |l| 
and a calculation that gives C5 ~ 96 to estimate an upper bound on the residual shape 
parameter, ri < 10"^ in lattice units. 



2. 3-hody untrapped ground state energy 

As a nontrivial test of the precision of our lattice method we have computed the lowest 
energy of three unitary fermions in a zero total momentum eigenstate; the energies of this 
state and higher eigenstates were computed to high accuracy by Pricoupenko and Castin 
in Ref. [29]. We performed the calculation for lattice sizes L = 8, 10, 12, 14, 16, tuning the 
coefficients of four operators for the L = 8 lattice, and five for the other lattices; for 
each lattice we used 1.5 — 1.9 x 10® scalar configurations. With a perfect one-body dispersion 
relation and this many two-body s-wave operators tuned, the leading L dependence of our 
result for the = 3 energy will be due to the untuned two-derivative two-body p-wave 
operator at 0{L~^)] subleading scaling would be due to the lowest dimension three-body 
operator, scaling as followed by the four derivative p-wave and d-wave two-body 

operators, scaling as L~^; for more details see 09]. In Fig. 4 we have plotted our results 
versus L^^ — the leading scaling behavior expected — including combined statistical and 
fitting systematic errors. Evidently the L > 10 numbers exhibit scaling nicely, while 
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FIG. 4: Energy of three untrapped unitary fermions in a zero total momentum eigenstate in units of 
the energy of three noninteracting fermions, E^ultrapped = 2 x (27r/L)^/(2M), plotted versus (bs/L)^ 
for L/bs = 8, 10, 12, 14, 16. The error bars include statistical and fitting systematic errors (for a 



discussion of these errors, see Sec. Ill A). The red band represents all possible two-parameter fits 



of the L/bs ^ 10 data to the function ci + C2/L^, refiecting both statistical and fitting systematic 
errors in our measurements, while the black line is the fit to the central values. The dashed line 
is the precise Pricoupenko-Castin result, Ref. ^29j, with which we agree to within our ~ 0.3% 
uncertainty. 

the L = 8 result is off, suggesting that L = 8 is too small a lattice to see the asymptotic 
scaling behavior. The red lines in Fig. 4 give the range of two-parameter fits of the L > 10 
data to Ci -|- C2/L^, which reflects the uncertainty in our data, while the black line is the 
fit of the central values of the data using the same fit function. At L — )■ cxd, the energy we 
obtain is 0.3735^°;°°^^ in units of the energy of three noninteracting fermions. As a result 
we find that our lattice action reproduces the Pricoupenko-Castin result to within our 0.3% 
uncertainty. 



D. External potentials 

Until now, we have concentrated on a system of interacting nonrelativistic fermions in the 
absence of an external potential. An external potential U may be introduced in a natural 
way by replacing the single particle interaction operator X defined in Eq. [5] with: 

X(r) ^ e-^/2^(r)e-^/2 ^ (29) 
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where the x matrix U is given by f/xx' = U{x.)Sx,x'- In the case of a harmonic trap, we 
use a potential of the form f/(x) = |kx^ centered about x = 0, and with simple harmonic 
oscillator (SHO) spring constant k. For fermions of mass M, the characteristic trap size is 



given by Lq = (kM)^^^^, and the oscillator frequency by w = a/ k/M. 

In the absence of interactions, the single fermion transfer matrix for our lattice theory is 
given by 

TsHo = e-pV^^^^e-^^^e-P'/^^^- , (30) 

which may be recognized as Trotter's product formula with 0{b^) time discretization er- 
rors'^. Specifically, temporal discretization errors are controlled by the dimensionless quan- 
tity (wfo,-)^, and are eliminated in the limit that — )■ in lattice units. 

Finite volume errors the other hand, are controlled by the dimensionless ratio L/Lq. In 
the continuum limit, finite volume errors for the noninteracting system may be computed 
analytically, since the SHO potential is separable. A plot of the energy dependence of 
the SHO on L/Lq is shown in Fig. |5]for several low energy single fermion states; at large 
L/Lq, the energies in units of u are just an integer plus the zero point energy 3/2 for a 
three-dimensional SHO. However, for very small volumes, the harmonic potential plays no 
role and the system is effectively a free particle in a finite box, with energies increasing 
proportional to ^ (^~qz^^ with decreasing L/Lq. The dashed lines in Fig. 5 indicate this 
limiting behavior for several SHO states. 

When tuned interactions are turned on, both temporal and spatial discretization errors 
are controlled in part by how the couplings are chosen. As was demonstrated in the previous 
section, by tuning the couplings one may completely eliminate both sources of discretization 
errors in the low end of the spectrum for two unitary fermions. Writing the transfer matrix 
for untrapped unitary fermions as 

r.ntrapped = 6"^^^^/^^'^ (1 - 6.V)e-''^P^/^^ = C''^^ (31) 

where "H is assumed to have been tuned free of discretization errors, then the transfer matrix 
for unitary fermions in a harmonic trap is given by 



The relation 7~(— tir) = T ^(^r) ensures that the energy can only suffer from corrections even in 6t-. 
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FIG. 5: L/Lq dependence of the SHO energies E^i 



nj ) corresponding to the single fermion 



states n = (0, 0, 0), (1, 0, 0), (1, 1, 0), (1, 1, 1), (2, 1, 0), (3, 1, 1), and (4,2,0). Sohd hnes indicate an 
exact continuum Hmit calculation, whereas the data-points indicate simulation results for uo = 0.005 
and Lq > 2. Dashed lines correspond to free fermions in a finite box (small L/Lq limit). 

^ g-6rC//2-fe?[C//2,p2/4M]+0(63)g-fe,-Hg-6rC//2+6?[C//2,pV4M]+0(63) 

^ ^-br('H+U)+Oibi) ^32) 

where {T-i + U) is the target Hamiltonian for trapped unitary fermions. We see that in the 
lattice definition of the trapped lattice Hamiltonian, "Htrappod = — ^ InTtrapped, temporal 
discretization errors appear at 0{h1). 

As was the case for noninteracting fermions in a harmonic trap, interacting fermions 
will possess spatial discretization and finite volume errors that scale as Bs/Lq and L/Lq, 
respectively. These errors must be explored numerically, and will be presented in detail in 
Sec. [ml 

III. ANALYSIS AND RESULTS 

In this section, we report results for the ground state energies of up to = 70 unitary 
fermions confined to a harmonic potential. We benchmark our method and systematic errors 
for up to = 6 against high precision solutions to the many-body Schrodinger equation, 
achieving agreement at 1%. We believe this is the first microscopic study to explore A^ > 6 
fermions in a trap without invoking a variational principle or requiring costly importance 
sampling. 

Numerical simulations of the trapped unitary Fermi gas have been performed with two 
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objectives in mind: evaluation of systematic errors using known few-body (A^ < 6) results as 
a benchmark, and numerical calculation of ground state energies of the many-body system 
(A^ < 70). We explore the question of whether one can use the trapped fermion data to 
extract the Bertsch parameter, defined as ^ = -£^untrappcd/-£^untrapped) where e';^]^^^^^^^ is the 
energy of noninteracting, untrapped fermions. ^ is related to the ground state energies of 
trapped fermions via the local density approximation [36] 

0{N-^^'')] (33) 



^trapped = 4tpedV^ (l ' ^^i^' " ^C^) (SiV)-^/^ 



where Ci and C2 are unknown phenomenological constants and -^trapped energy of 

noninteracting trapped fermions. 



trapped " ^ ^ ■ \AV 

Note that if (ci — |c2) — 1, then for = 70 and ^ ^ 0.4 one finds the subleading term 
in the expansion Eq. ( [33| ) to be the same size as the leading term, suggesting that A^ = 70 
is not enough particles for the trapped system to be considered near the thermodynamic 
limit. In fact, that is what we find: we see significant shell structure all the way up to 
A^ = 70 and conclude that we are not yet in the thermodynamic limit. This is in contrast 
with what we find in the untrapped case, where shell structure disappears at much lower 
A^ |l9]. At A^ = 70 the system has not yet reached the thermodynamic limit, we are not 
able to extract the value of ^ or the unknown parameters ci and C2 from the trapped data. 
Our data does, however, give information about possible differences in how the trapped and 
untrapped systems approach the thermodynamic limit. 



A. Extraction of ground state energies 

The energies of multi-fermion systems may be extracted from correlation functions using 
conventional techniques. Given a correlator C(r) describing the Euclidean time evolution of 
some N-fermion initial state (source) at time slice zero into some final state (sink) at time 
slice r, a generalized effective mass may be defined as 



At 



log 



C(r) 



C(r + Ar) 



(35) 



which satisfies limr^oo ^eff{T) = Eq, where Eq is the ground state energy of the system. 
At late times, energies are given by a plateau in the effective mass, with excited state 
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contamination falling off exponentially in the energy difference between lowest and first 
excited states. For noisy correlators, a stride of Ar > 1 may be used to facilitate detection 
of the time window over which a plateau appears. 

For large numbers of fermions, the standard effective mass exhibits a distribution overlap 
problem (see Appendix [C]). For this reason, we utilize the effective mass defined using the 
cumulant expansion truncated at C(A^k), 




n=l 



where Kn(r) is the nth cumulant of log(C(r)). Details of this technique may be found in 
Appendix [Cj and details of the particular strategy used for systems of trapped fermions 
will be discussed in Sec. IIII D 1[ 

To extract the energies of the system, we perform correlated fif s to the plateau region 
of the effective mass associated with the fermion correlator. Statistical error estimates 
are obtained by resampling the data using the bootstrapping technique. Fitting systematic 
errors are found by varying the endpoints of the fitting interval. For small A^ ^ 8, contami- 
nation from excited states persists to very large Euclidean times. Because of this, the data 
we fit is quite noisy and determining the plateau region becomes difficult. For this reason, 
we vary the endpoints of our fits by 6t = ±10 to account for any systematic error due to 
the choice of fitting region. For large A^ we find that it is sufficient to vary the endpoints 
of the fit region by an amount 6t = ±2 to determine our fitting systematic errors. Because 
the plateaus are well-resolved for many time steps we do not find significant deviations in 
the error bars by considering larger variations of the endpoints. 

B. Ensembles and parameters 

A complete analysis of the systematic errors due to finite volume and lattice spacing 
artifacts requires performing scans in the parameters L, Lq and u. Since performing such 
scans would be prohibitively costly for large numbers of fermions, we have instead chosen 
to generate two sets of ensembles that allow us to address these questions in a cost-effective 
manner. 

The first set of ensembles consists of a series of scans in the aforementioned parameters, 
while restricting the number of particles to values A^ < 6. Restricting the number of fermions 
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TABLE III: Many-fermion simulation parameters for trapped fermion using the pairing wave func- 
tion given by Eq. B6 with /3 = l/(\/2Lo). 
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in this way greatly reduces the computational resources required, and also permits a higher 
resolution in the parameter scans. Few fermion ensembles of size N^onf = IM were generated 
for L = 48 and L = 64 lattices using trap sizes Lq = 3, 4, 5, 6, 6.5, 7, 7.5 and 8. Scans were 
primarily performed at = 0.005, however, several simulations were also performed at 
u = 0.01. The temporal extent for all of the few-body lattices was T = 80. 

With guidance from our analysis of the systematic errors of the few-body system, we then 
performed a more targeted set of simulations for up to = 70 fermions, using parameter 
choices L = 48, 54 and 64, Lq = 7 and 8, and u = 0.005. The parameter choices used in our 



< 70 simulations are detailed in Table III For our simulations of up to = 70 fermions, 
we have generated approximately one million configurations for each value of the volume 
and trap size, using a total of less than one million CPU hours. In all of the trapped fermion 
studies, we have used No = 4 tuned couplings for the interaction. 

Details of our construction of multi-fermion correlation functions are given in Appendix [B| 



Following |32], we use a modified Slater determinant Eq. B3 and Eq. B4 to include pairing 
correlations. The sinks are constructed from the two-particle wave functions defined in 



Eq. |B6[ For all correlation functions, the free parameter appearing in Eq. B6 was chosen as 
(3 = I/v^Lq. Multi-fermion sources where constructed from free SHO single particle wave 



functions |n°") with a = (4)t) provided in Table VI, with i < N/2. The sources involving 
odd A^ for our few fermion studies were obtained by removing a single fermion from the 
highest shell, as described at the end of Sec. [Bj 
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C. Few-body Results 



To reach the continuum and infinite volume hmits we require bg Lq <^ L, and b^- ^ 
1/u. To balance the need for small temporal discretization errors with the computational 
cost associated with the number of time steps required to reach the ground state, we have 
chosen ubr = 0.005 for this study. For small A^, we find that the discrepancies in the energies 
for ubr in the range 0.005 — 0.01 are about 0.5%, and are within our error bars. 

For a given box size, the choice of Lq must take into account both discretization errors and 
finite volume errors. The expectation is that for small L^/bg spatial discretization errors will 
dominate. The discretization is implemented as a hard cutoff in momentum space, which 
may be interpreted as an infinite potential at the edge of the Brillouin zone. Sensitivity of 
the state to this infinite potential results in an increase in the associated energy. Conversely, 
for large Lo/bg and fixed volume, finite volume errors will dominate. The periodic boundary 
conditions in space result in attractive interactions from image particles, causing a decrease 
in energy. Thus, measurement of the ground state energies as a function of L^/bg and L/Lq 
is necessary to determine at which value we can minimize both types of error. 

Fig. [6] presents our findings for the ground state energies of = 3, 4, 5 and 6 fermions, 
with Lo/bg ranging from 3 — 8 and fixed L/bg=A8. Also indicated in this figure are the 
ground state energies for unitary fermions in a trap quoted in |57] , which were obtained by 
numerically solving the multi-fermion Schrodinger equation using the correlated gaussian 
(CG) method. Using the results of [ST] as a benchmark, we find that for L^/bg < 7.0 our 
discretization errors are significant. Above this value, however, we find that the extracted 
energies are independent of Lo/bg, indicating negligible discretization errors in this regime^. 

In Fig. [7| we present the L/Lq dependence of the energies for A^ = 3,4,5 and 6, with 
L/bg = 48 and 64 and Lq > 7. The consistency of the results between the different volumes 
indicates that finite volume errors are negligible within statistical uncertainties for L/Lq 
ranging from 6 — 9. The good agreement of our Lq > 7 data in Fig. [7] with the benchmark 

^ In [35], we found that our results agreed with those of [57] for values of Lq sa 4. However, it became 
evident that this agreement resulted from a delicate cancellation between temporal and finite volume 
errors, and that each source of error was individually rather significant. In this work, we have reduced 
the temporal discretization errors with an improved form of the potential; this improvement results in 
temporal errors appearing at an order higher in ubr- We have also chosen a smaller value for ubr, and 
checked that the results are consistent for both the smaller (ojbr = 0.005) and larger {wbr = 0.01) values. 
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FIG. 6: Ground state energies (in units of uj) as a function of L^/bg at fixed L/6s = 48 for various 
values of N . Dashed lines are results from [57]. 



TABLE IV: Results for ii'trappcd/'^ for ^ 6, including combined statistical and fitting systematic 
errors (first row). For comparison we give the exact = 3 result [15j and results of Ref. |57) 
(second and third rows). 
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4.273(2) 
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energy values indicates the absence of any residual errors. We have performed a constant 
correlated fit using all the data in Fig. [7] to obtain infinite volume, vanishing lattice spacing 
results for the few particle energies. Our final fit results for < 6 are indicated in Fig. [7[ 
and presented along with a comparison to an exact result for A^ = 3 [I5] and high-precision 



Hamiltonian results of [57] for A^ = 4 — 6 in Table IV 



For A^ > 6, it is likely that both discretization and finite volume errors will grow, since we 
expect the wave function to spread out in both position and momentum space when more 
particles are added to the system. Numerical evidence suggests that extrapolations become 
necessary for A^ > 20. More details of our analysis for larger A^ will be presented in the next 
section. 
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FIG. 7: Fit results for the ground state energies (in units of cj) for < 6 as a function of L/Lq 
for two volumes, L = 48, 64, and three trap sizes, Lq = 7.0, 7.5, 8.0. The result of a correlated fit 
in N, T of all data is shown as a red line, with a red band showing the combined statistical and 
systematic errors. The results from [15] (N=3) and [57] (N=4-6) are given by dashed lines, with 
any associated error bars shown by hatched regions. 

D. Many-body Results 

1. Statistics 



Examples of an effective mass plot obtained using the conventional definition, Eq. |35 
with At = 1, for = 30 and = 70, are shown in Fig. |8| Although we have found good 
signals for most values of A^ at short times, at later times the effective mass plot shows clear 
evidence of a distribution overlap problem. 

Generally speaking, since the sources and sinks used to compute these correlation func- 
tions are different (see Appendix [B]), positivity of the correlator is not guaranteed. Further- 
more, there is no reason to expect that effective masses obtained from them will decrease 
monotonically as a function of time. When analyzing effective mass plots, one must therefore 
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FIG. 8: Effective mass as a function of r (lattice units) for = 30 and = 70 (L = 48, Lq = 8.0) 



using the conventional definition for the effective mass Eq. 35, Fits are represented by horizontal 
bands. 

be capable of distinguishing between local minima in the effective mass and true plateaus in 
order to extract reliable ground state energies. In cases like Fig. [s] (iV = 30) where a plateau 
begins (r ~ 6) well before the onset of a distribution overlap problem (r ~ 22), one may 
easily make the distinction between local minima and true ground state plateaus. However, 
in cases like Fig. [s] (A^ = 70) where the onset of an overlap problem and the beginning of a 
(potential) plateau coincide, the distinction becomes less clear. 

In light of the considerations above, effective mass calculations based on the cumulant 



expansion, Eq. 36 (for details, see Appendix [C]), were used to help establish and extend 
plateaus into the region of poor overlap. The result of this technique is demonstrated in 
Fig. |9] for = 70 with ranging from 3 — 5. Comparing the effective mass obtained 
from the cumulant expansion with that obtained by conventional means, we see a marked 
reduction in the overlap problem for r > 12, although the noise at late times increases with 
N,^. Nonetheless, the extension of the plateau region to larger times gives us confidence that 
we have reached the ground state, and allows us to perform fits over much longer temporal 
extents. 

Examples of the fits obtained using the cumulant expansion for small (A^ = 10), moderate 
(A^ = 30), and large (A^ = 70) numbers of fermions are shown in Fig. [To] For small A^, 
the plateaus appear at late times where we find that the cumulant expansion converges 
slowly. However, because the overlap problem is less severe for small A^ versus large A^, we 
may corroborate our cumulant results with those obtained using the conventional effective 
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FIG. 9: Effective mass as a function of r (lattice units) for N = 70 (L = 48, Lq = 8.0) using 



the cumulant expansion (Eq. 36) for up to N,^ cumulants. The fits from each are represented by 
horizontal bands. The conventional effective mass is shown in gray. 



mass (red bands). For = 10, convergence occurs at N,^ ^ 6, whereas for all > 12, 
convergence occurs at A^^ ~ 3. Higher cumulants may be used to further extend the plateau 
region without significant growth in the error bars for the fits. This is presumably because 
the fit results are highly influenced by the early time portion of the plateau region, where 
the errors on the effective mass are relatively unaffected by an increase in the number of 
cumulants. The nearly exact agreement between results obtained using = 3,4, 5 leads us 
to conclude that systematic errors introduced by truncation of the cumulant expansion are 
negligible. 



2. Systematics 

To account for systematic errors arising from finite volume and lattice spacing effects, 
we have performed the calculation for three volumes {L = 48, 54, 64) and at two values 
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FIG. 10: Fit results for the ground state energies (in units of uj) for small, moderate, and large A'^ 



(L = 54, Lq = 8.0) using the cumulant expansion including A^^^ cumulants (Eq. 36). The ground 
state energy extracted from the conventional effective mass is given as a red band. 

of the trap size (Lq = 7.5,8.0). We find that as more particles are added to the system, 
the discrepancies between results at different volumes grows. The dependence on the lattice 
spacing is less clear, particularly due to the fact that changing Lq changes not only the lattice 
spacing dependence {hg/ Lq), but also the finite volume dependence {Lq/ L). To separate these 
effects, an infinite volume extrapolation was performed for each value of Lq using correlated 
fits of the data to 



(37) 



over the plateau regions of each effective mass plot. This form of the extrapolation function 
is a simplified version of the ansatz that finite volume errors depend on the probability, 

/oo 
\ij{x/LQ)\Hx (38) 

that the ground state wavefunction extends outside the box. We also make use of the fact 
that for unitary fermions, iP{x/Lq) is given asymptotically by a direct product of noninter- 
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FIG. 11: Volume dependence of the ground state energies (in units of uj) for moderate {N = 24). 
The data points indicate the individual results for our six values of L/Lq. An infinite volume 
extrapolation is shown as a solid line, while a band represents the associated statistical and fitting 
systematic error bars of the extrapolation. The upper plots show separate fits to the Lq = 7.5 and 
Lq = 8.0 points using Eq. 37 The lower left plot shows a combined fit using both values of Lq, 



and the lower right plot shows the combined fit with a final error band obtained by combining the 
statistical and fitting systematic errors from all three extrapolations. 

acting harmonic oscillator wavefunctions. We find that including wavefunctions from higher 
shells results in negligible change from the infinite volume extrapolations obtained using 
Gaussian fits. These differences may ultimately be absorbed into the fitting coefficients 
A,B. 

For the two different values of Lq we find that the infinite volume extrapolations are 
consistent within error bars, indicating that spatial discretization errors are smaller than 
the combined statistical, fitting systematic, and infinite volume extrapolation errors. For 



this reason, a third fit was also performed to all six data points simultaneously (see Figs. 11 



12). The spread between the three fits gives an approximation for any remaining spatial 
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FIG. 12: Volume dependence of the ground state energies (in units of w) for large (N = 70). 
The data points indicate the individual results for our six values of L/Lq. An infinite volume 
extrapolation is shown as a solid line, while a band represents the associated statistical and fitting 
systematic error bars of the extrapolation. The upper plots show separate fits to the Lq = 7.5 and 
Lq = 8.0 points using Eq. 37 The lower left plot shows a combined fit using both values of Lq, 



and the lower right plot shows the combined fit with a final error band obtained by combining the 
statistical and fitting systematic errors from all three extrapolations. 

discretization errors. For the final result, we added the statistical and fitting systematic 
errors from each fit in quadrature individually, and used the outer envelope to represent our 
total statistical, fitting systematic, extrapolation, and spatial discretization error. 

3. Final Results 



Our results for the energies in units of u and their corresponding errors are reported 



in Table |V| In Fig. 13 we plot the results for the ground state energies in units of the 
energies for the corresponding noninteracting system, -EtrappedZ-E^^^ppg^j. For comparison, we 
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TABLE V: Ground state energies as a function of N in units of lj. The error represents the combined 



statistical, fitting systematic, finite volume, and spatial discretization errors. See Sec. HID 4 for 
possible additional systematic errors. 
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also show the results from two fixed-node calculations: a Green's function Monte Carlo 
(GFMC) approach [12] and a Diffusion Monte Carlo (FN-DMC) approach |41j. By using 
the fixed-node constraint along with a variational principle, both of these methods provide 
upper bounds on the ground state energies. We find that our energies are consistently lower 
than those obtained using both of these methods. Interestingly, fixed node calculations do 
not display the shell structure which is clearly present in our data. It is evident that this 
shell structure diminishes for large A^, where eventually the thermodynamic limit should be 
reached. 
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FIG. 13: Ground-state energies of trapped unitary fermions in units of the corresponding energies 
of trapped noninteracting fermions as a function of A^. For comparison, we show results from 
GFMC US], FN-DMC glj, and CG ^ methods. 

4- Possible additional sources of systematic error 

To calculate the error bars quoted in Table |V] we have taken into account statistical, 
fitting systematic, extrapolation, and lattice errors. We note additionally that the spacing 
between the energy levels associated with breathing modes [18], 2ubr = 0.010, is smaller 
than the inverse temporal extent of our lattice (1/T ^ 0.017), but larger than our quoted 
error bars. Furthermore, as an increasing number of particles are added to the system, a 
near continuum of different angular momentum states may result, also of 0{u)br). 

These excited state contributions could lead to systematic effects due to a failure to reach 
the ground state of the system. If excited state contamination is present in our results, it 
is possible that the overlap of our chosen sources and sinks (see Appendix [B]) with these 
excited states is shell dependent, causing our results to exhibit shell dependence even if 



this is not a property of the ground state. However, as noted in the beginning of Sec. Ill 
we do not have any reason to believe we are near the thermodynamic limit, so it is quite 
conceivable that the shell structure we observe is a physical property of the ground state. 

The energy splittings are the same size for small as for large A^, thus we might expect 
that if we are able to see the ground state within our time extent for small A^, the same 
could be true for large A^. Because our results for small A^ agree with those from benchmark 
calculations, we can be assured that we have found the ground state in this case. This 
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FIG. 14: Effective mass as a function of r (lattice units) using two different sources for a half-filled 
shell (A'^ = 14) and a nearly closed shell {N = 18) within the second shell (L = 48, Lq = 8.0). The 



effective mass was calculated using the cumulant expansion with A^^ = 6 using Eq. 36 The blue 



circles were generated using a source constructed by filling the single-particle states in the order 



given in Table VI while the red stars were generated using a source constructed of random linear 
combinations of the single-particle states within each shell. 
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FIG. 15: Effective mass as a function of r (lattice units) using two different sources for a closed 
sheU {N = 20) and a half-ffhed shell {N = 30) within the third sheU (L = 48, Lq = 8.0). The 
effective mass was calculated using the cumulant expansion up to A^^^ = 6 using Eq. |36} The blue 
circles were generated using a source constructed by filling the single-particle states in the order 



given in Table VI while the red stars were generated using a source constructed of random linear 



combinations of the single-particle states within each shell. 

implies that the wavefunction overlap with excited states is very small. 

To better quantify any possible effects from excited states, we may consider a correlation 



33 





74 




72 




70 


O) 




a 
a 


68 


b 










66 




64 




67. 



i 


N=30 


* 








= Early time fit (random source) 




— Early time fit (pure source) 




* Late time fit (random source) 




• Late time fit (pure source) 



FIG. 16: Fit results for the ground state energy using the two sources shown in Fig. 15 (right). 
The bands show results from fitting each source at early times (r ~ 5), while the data points show 
fit results as a function of the number of cumulants included for late times (r ~ 30). The late 
time fits from both sources are approaching the early time fit for the pure source, while the early 
and late time fits for the random source do not agree, indicating that the random source does not 
reach the ground state until later times. 

function whose long time behavior is dominated by two terms, the first corresponding to the 
ground state, the second to a breathing mode 



{Eo+2uj)t 



(39) 



where Zq and Zi represent the overlaps between our sources and sinks with the true ground 
and breathing mode states, respectively. Recall that the signs of Zq and Zi need not be 
positive due to the use of inequivalent sources and sinks. 

For large we typically find a plateau for time ranges r ~ 5 — 30. If we assume equal 
coupling of our sources to the ground state and breathing mode {Zq = Zi), this would 
contribute to a drift in the effective mass plot of about O.la; for the time range considered. 
This is of approximately the same size as our statistical error bars in this region, so it is 
conceivable that such a drift would not be detected. 

One possible test to detect contamination from excited states is to repeat the calculation 
using a source that consists of random linear combinations of the states within each shell 



from Table |VI[ Using Eq. [39| one may show that the shift in the ground state energy for 
small UT is approximately 



Eo + 2uj 



34 



(40) 



If there is contamination from the second term and the new source changes the overlap with 
the excited state by at least 0{1), this would give a shift in the effective mass plot of 0{u) 
for all times considered (r < 64). 

We find that the effective mass plots produced using the random sources agree with those 



for our original ("pure") sources for within the first two shells (see Fig. 14). For the third 



shell (Fig. [15| ), the effective mass for the random source begins at higher values for both 
= 20 (closed shell) and = 30 (half-filled shell), however, the two sources begin to agree 
around r ~ 30. 

By performing fits using the cumulant method, we find that in fact the random source 
plateaus at a later time than the pure source; results from fitting both sources at late 



times (r ~ 30) are consistent with each other (see Fig. 16). While the cumulant expansion 
converges too slowly at late times for us to extract a reliable ground state from these fits, it 
is clear that the results from both sources are approaching the early time (r ~ 5) fit for the 
pure source, giving us confidence in the energies extracted from this source. 

Thus, the random source test supports a lack of contamination from excited states in our 
quoted results. However, we do note that there is no guarantee that randomizing the source 
changes the Z-factors by at least 0{1). Further analysis will be necessary to definitively 
estabhsh this conjecture. 



IV. CONCLUSION 



We have developed a new lattice method for studying large numbers of fermions at 
unitarity. The action is highly improved, so that our results require no extrapolation to zero 
range. In addition, we've applied a new method for calculating correlators from long-tailed 
distributions |1S], through which we are able to evade costly importance sampling. Our 
results agree with those from high precision solutions to the Schrodinger equation for < 6 
trapped fermions [57j, as well as with the energy of A^ = 3 untrapped fermions calculated 
by Pricoupenko and Castin [29]. Due to the low cost of the simulation we are able to then 
study up to A^ = 70 trapped fermions, finding lower values than pubhshed results. One 
feature we find is that shell effects persist at the ~ 2% above A^ = 40 fermions, making it 
impossible to extract a reliable value for the Bertsch parameter ^. The shell effects we find 
are much more pronounced than what we see for untrapped fermions ^9] . 
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In a future work we will present results for the homogeneous system of up to N — 66 
unitary fermions in a box, including our extraction of the Bertsch parameter, ^, as well as 
data on the superfiuid gap and integrated contact density for unitary fermions in a box. 

We believe this method could be applicable for a wide variety of nonrelativistic many-body 
systems, and these studies of unitary fermions, in addition to their inherent value, pave the 
way for investigations of more complex systems at zero temperature. 
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Appendix A: Tuning 

The two particle transfer matrix is a linear function of the two-body couplings C2n' 



where 7/ree is the free fermion transfer matrix, and 72n contains contributions to the in- 
teraction C starting at order 2n in momenta. Eigenvalues Afc(C) of this transfer matrix, 
however, are nonlinear functions of the couplings. We may compute the derivative of these 
eigenvalues with respect to the couplings using the Feynman-Hellman theorem: 



No-l 



(Al) 



n=0 



dXk 



{'ipk\T2n\-ipk) , 



(A2) 



dC2n 
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where T(C)|V'fc) = Afe(C)|V'fc) for k = 1, . . . , dim(T) and the eigenstates \ipk) imphcitly 
depend on C. Tuned values for the Nq couphngs are defined as the values C2„ {n = 
0, . . . , No — 1) for which x^(C) = 0, where 

No 

x\C) = J^5\,iCf (A3) 
fe=i 

and 5Afc(C) = Ajt(C)/A^, — 1. Starting from an initial guess for the couplings C2„, we may 
iteratively search for the solution to = using: 



provided the inverse of the A^^cj-dimensional sub-matrix W exists, where Wkn = for 
n = 0, . . . , Nq — 1 and k = 1, . . . , Nq. The small parameter e may be chosen adaptively in 
order to improve the convergence of the iterative procedure. 

Appendix B: Correlation Functions 

Multi-fermion sources may be constructed from direct products of single particle states 
where i = 1, . . . ,N„ labels each state with quantum number a and a = (t^i) labels 
the species. In order to satisfy Fermi-Dirac statistics, fermions of the same species must 
have different quantum numbers. As is well-known from quantum mechanics, a simple way 
to impose the proper anti-symmetrization requirements on multi-fermion states is to use 
Slater-determinants. Thus correlation functions of N = N^ + N^ fermions may be expressed 
as: 



where 5"^ is an Ai'o-- dimensional Slater matrix corresponding to the species a, given by 



Although it is not a requirement, a convenient choice for the single particle states jaf) is 
to use eigenstates of the noninteracting system. For trapped fermions, they are SHO states 
{a = n) in the Cartesian basis. A list of the sources used in our simulations is provided in 




(A4) 




(Bl) 




(B2) 



Table |VD 
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TABLE VI: Single fermion sources used in trapped (a = n) fermion calculations. 
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Typically multi-particle sources constructed from single particle states possess poor over- 
lap with the unitary Fermi gas ground state. This may easily be seen from the fact that 
at early times, where few interactions have occurred, the correlation function falls off ex- 
ponentially like that of free fermions with a Z-factor near unity. A better approach is to 
incorporate pairing correlations into the interpolating field by constructing sources and sinks 
out of two-fermion wave functions ^32j. In practice, such an approach may only be carried 
out at the sink, however, because our numerical approach requires that sources be separable 
functions; this is none-the-less adequate to achieve far superior overlap with the ground 
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state. A consequence of using sources and sinks that differ is that correlation functions and 
effective masses need not be monotically decreasing functions of time. Thus when studying 
correlators of this form, care must be taken to distinguish shallow local minima in effective 
masses from true plateaus. 

For A^''^ = N-^ = N/2, these considerations lead us to study correlation functions of the 
form: 

CN„N,{r) = {detS^\T)) , (B3) 

where 

Sj]{T) = {m-\r,0)^K-\r,0)\aja]) (B4) 

and |Q;~''a^) = \a^) (g) \a^)- In the coordinate basis, we consider two-fermion states |\t') of 
the form (x~''x''^|\E') = \E'(r,.e;) where Trei = x-l- — x''^ is the relative coordinate of the two 
fermions. It is helpful to express the two-particle wave functions as a Fourier transform: 
\E'(rre«) = / (ip^(p)e~P'''''=', allowing Eq. 



B4 



to be written as 



Sf]i{T) = J2nP){p\K-\r,0)\ai){-p\K-\TMc^]) . (B5) 
p 

Since the projection onto the sink involves only a single sum over momenta, evaluation of 



Eq. B4 scales like 0{L^) rather than the usual 0{L^), 

Numerical evidence suggests that the best choice for \E'(rrez) is a lattice approximation to 
the two-particle s-wave solution to the continuum Schrodinger equation for unitary fermions, 
which possess a l/lr^ezl singularity. We therefore consider a momentum space wave-function 
of the form 

^iP) = ^d(^-^] , (B6) 



IpI V2/3, 

where d{x) is Dawson's integral function. Note that the wave-function has a free parameter 
P which may be tuned to maximize the overlap with the ground state. Physically one expects 
/3 ~ 1/a/2Lo, and this is what we use in practice. 

For odd numbers of fermions, such as in our few-body studies, one may construct a mixed 
matrix built out of both single and two fermion wave functions. In the case A^^'' = A^^ + 1, 
one may may construct such a Slater matrix by replacing row i of 5*^''^ with the same row of 
S'^. This replacement corresponds to a removal of the i-th single fermion state aj from the 
source and thus also breaking the pair involving state i at the sink. 
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FIG. 17: Natural logarithm of the N=4 fermion correlation function for untrapped unitary fermions 
of mass M = 5 on an L = 10 lattice as a function of sample size. Dashed line indicates the result 
obtained using N^onf = 2S configurations. 

Appendix C: Measurement Strategy 

Our studies have shown that for large numbers of fermions, the effective mass, 

C(r) 



(CI) 



_C{t + At)\ ' 

obtained from correlators measured in the conventional way, as a sample average of propaga- 
tors measured on background (p configurations, often exhibits a distribution overlap problem. 



Several indicators for this problem include: 1) lack of 1/ ^jNconf scaling in the error bars for 
correlators and measured quantities derived from them, 2) sudden jumps in estimated quan- 
tities and their associated error bars as a function of the sample size N^onf and 3) continued 



growth of effective masses at late times, with no evidence for a plateau. Fig. 17 provides a 
mild demonstration of the third case for = 4 untrapped unitary fermions; plotted is the 
logarithm of the correlation function C(r) which has been estimated using ensembles of size 
ranging from Nconf = 0.06M to 3.84M configurations. An estimate of the effective mass at 
late times, quantified by the slope of — logC(t) in this figure, appears to decrease as the size 
of the ensemble is increased. As the number of configurations in the ensemble is increased 
by several orders of magnitude, C(r) eventually appears to converge to what is expected to 
be the true value of the correlator, indicated by the dashed line and estimated using a much 
larger ensemble of size Nconf = 2-B configurations. 

In order to understand this behavior better, it is instructive to study the distribution of 
the correlator as a function of r. The distribution of an arbitrary operator Y[(j)) measured 
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FIG. 18: = 4 fermion correlator and natural log-correlator distributions at various time separa- 
tions r for unitary fermions of mass M = 5 on an L = 10 lattice. Solid curves in the log-correlator 
distribution plot correspond to Gaussian fits to the distribution. 



on a background field configuration is given by: 

p{y)= l[dct>]p{my{<p)-y) ■ 



(C2) 



A plot of the correlator distribution, taking Y[(j)) = C^{t) and y = c, is shown in Fig. 18 for 
= 4 fermions at several values of r, and demonstrates the formation of a long tail in the 
late time limit. Also shown is a corresponding plot of the distribution for the logarithm of 
the correlation function, taking Y{(f)) = logC0(r) and y = logc, along with the results from 
a Gaussian fit to the histograms. The excellent agreement between the fit results and the 
measured log-correlator distribution suggests that the multi-fermion correlation function is 
log-normally distributed, or nearly so. Such distributions are known to possess very long 
tails which dominate the distribution mean, and undersampling the tail can result in an 
underestimate in the correlation function, and thus an overestimate in the energy obtained 



from Eq. CI at large times, as was evident in Fig. 17 

Provided we know the underlying distribution for the correlation function, we may esti- 
mate the number of configurations required such that the sample average C(r) is normally 
distributed. Deviations of the sample mean from the normal distribution may result in an 
overlap problem and reflect the fact that the sample size is too small for the central limit 
theorem to apply. In particular, if a and p are the second and third central moments of 
the correlator distribution function, then by the Berry-Esseen theorem, one should show 
that the condition Nconf ^ P^/'^^ holds before invoking the central limit theorem. This 
condition comes from quantifying the deviation in the cumulative distribution function for 
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FIG. 19: Plot of the cumulative distribution function for (C — (C)) Nconf I ^ for C estimated using 
ensembles each of size Nconf = lOOX; (C) is estimated using 2B configuration. The dashed 
line is the cumulative distribution function for the standard normal distribution. 



{C — (C)) a/ Nconf I o" from that of the standard normal distribution, where C is the sample 
mean of the correlation function obtained from a sample size Nconf- An example of the cu- 
mulative distribution function of this quantity for = 4 unitary fermions at several times is 
shown in Fig. 19, and was obtained from 20K ensembles each of sample size Nconf = lOO-ft'. 
The true mean (C) was estimated using a sample of size Nconf = 2-B configurations. In 
this example, we find that for r = 24 there is little deviation from the standard normal 
cumulative distribution function, whereas for r = 36, significant deviation is evident. 

In [IH] , it was shown within mean field theory that the log-correlator distribution function 



defined by Eq. C2 is Gaussian with mean y = log Z + Eq{N)t and variance = ^Eo{N)t, 
where Eq{N) is the free gas ground state energy for noninteracting fermions {N/2 fermions 
of each species) and log Z is the corresponding overlap between the ground state wave func- 
tion and source and sink wave functions. This in turn implies that the correlator distribution 
is log- normally distributed in mean field limit. We may use the Berry-Esseen theorem along 
with our mean field result for the correlator distribution to estimate the minimal number 
of configurations required for a given value of and r. The results is Nconf ^ e^^f ^o(^)'^, 
scaling exponentially in the time and free gas energy. Applying this result to the case = 4, 
we one finds that Nconf ^ Si^' configurations are required for r = 24, Nconf ^ 175K for 
r = 36, and Nconf ^ lOM for r = 48. The onset of an overlap problem around r ~ 32 in 
Fig. 19 obtained from Nconf = lOOK configurations is consistent with the prediction based 



on the application of Berry-Esseen theorem to our mean field calculation. 
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The traditional technique for avoiding difficulties associated with distribution overlap 
problems is to use importance sampling in the Monte Carlo simulation. In the case of large 
numbers of fermions, this might be achieved by reweighting the probability measure by either 
the correlation function at some late time, or some other carefully chosen weight factor. In 
the former case, one might use the product p(0)C0(ro) for an arbitrary but large value of tq as 
a probability measure for the auxiliary fields, and then measure ensemble averages of the ratio 
C0(r)/C^(To) to estimate the correlator at times r. In taking such an approach, however, 
the ensembles generated are typically only suitable for estimating a specfic operator (e.g., a 
single correlator at a specific value of A^) or a small class of operators, and are inappropriate 
for most others. Consequently, the simulation cost is enhanced by the number of operators 
being measured in addition to the difficulty of performing unquenched simulations using a 
far more complicated effective action for the auxiliary field, which generally will involve the 
logarithm of a correlation function. This may be likened to performing a simulation in the 
Grand Canonical ensemble, where a different simulation must be performed at each value 
of chemical potential to achieve estimates of the energy as a function of density. 

A far more efficient approach proposed in |18] is to find a better estimator for C(r) 
that is free from the distribution overlap problem rather than rely on importance sampling. 
Provided C^{t) > for every 0, a systematic method for extracting useful information from 
an undersampled log-normal-like distribution may be devised by considering the cumulant 
expansion: 

logC^„.(r)^^^, (C3) 

n=l 

where K„(r) is the n-th cumulant of the distribution for log C^{t), which is presumed to be 
nearly normally distributed. In this expansion, systematic uncertainties associated with the 
truncation of the series at order are traded for statistical uncertainties associated with 
including increasing numbers of cumulants which have been estimated from an ensemble of 



finite size. For a perfect log-normally distributed C^{t), Eq. C3 is exact at A^^ = 2, since all 



Since effective masses depend only on tlie ratio C(t + 1)/C(r), tfie overall normalization of correlation 
functions determined from an ensemble average of 1/C^{to) using p{(j))C^{TQ) as a probability measure is 
unimportant. 



■^^ For the case N^^ — N^, one can show explicitly that correlators of the type defined in Eq. B2 and Eq. 
are positive for every background field configuration. 
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higher order cumulants vanish. In practice, if the correlator distribution is not log-normal, 



deviations in the distribution would be quantified by the non-zero contributions to Eq. C3 



from Kn with n > 2. Such contributions-one would hope-are relatively small, allowing one 
to reliably obtain an estimate for log C(r) based on estimates of k„. 



The generalized effective mass associated with each partial sum in Eq. |C3| may be ex- 
pressed as: 



AN.), 



mlT^'ir) = i- ^ 1 [ft:„(r) - k^t + At)] . (C4) 



n=l 



By studying Eq. C4 as a function of A^^, one may determine the ideal value A^* for which the 
statistical uncertainties and truncation errors become comparable. Such an A^* then defines 
a best estimate value for the effective mass at a given time r. Alternatively, we may define 
an energy E^^ = limT-_j.oo ^^^ff{'^Y'^ and study its convergence as a function of A^^. In all of 
our studies, we use the latter approach. 

Finally we comment on the applicability of the cumulant method to odd numbers of 
fermions. In the case A^^ = A^^'' + 1, analysis using the cumulant expansion breaks down, 
since negative correlators C<^(t) may exist. For large numbers of fermions, we find that 
the fraction of negative correlators in a given ensemble is typically less than a few percent, 
however. Furthermore, unlike the positive part of the distribution, the negative part exhibits 
no long tail at large r. This suggests that the positive and negative parts of the distribution 
may be treated not only independently, but also differently: for the positive portion one may 
use the cumulant expansion technique, and for the negative portion a standard ensemble 
average, and the results may then be combined. Although we do not consider odd numbers 
of fermions in this paper, we believe these considerations will be of importance in future 
studies of the pairing gap, which requires accurate estimates of the energy for both even and 
odd A^. 



Appendix D: Simulation details 

The lattice theory described in Sec. [TT] has been implemented on a number of clusters 
and massively parallel architectures. As a result of the low memory footprint, which scales 



Although we have not proved the convergence of m^jj- (t) as a function of t, aU of of our numerical 
evidence suggests that this quantity tends to a constant at late times. 
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like 0{N X L^), and extremely fast character of the algorithm |13], we use an embarrass- 
ingly parallel implementation: multiple streams are farmed out to many different cores and 
random generators associated with each core are seeded independently of each other. Sim- 
ulations were performed in double precision, and random numbers were generated using 
Liischer's Ranlux pseudo-random number generator [59] • Preliminary studies have shown 
no statistical advantage to using Gaussian auxiliary fields over Z2 noise, and therefore all of 
our studies have been performed using the latter. 

Due to the extremely fast nature of our algorithm, it was necessary to perform ensemble- 
averages of multi-fermion correlator data on-line in order to eliminate bottlenecks associated 
with file I/O and to also reduce storage requirements. Data was therefore ensemble averaged 
into Nb blocks of size Nconf/Njs, where Nconf is the total number of configurations generated. 

was chosen small enough to avoid the I/O and storage issues, but also large enough to 
maintain adequate control over the statistical errors in our analysis. 

We have checked our algorithm and implementation by comparing numerical predictions 
of the energies for several exactly soluble systems with their known solutions. The four- 
fermion interaction, for instance, was checked using a high precision measurement {Nconf = 
4.7 B configurations) of the lowest and first excited state energies for two unitary fermions 
of mass M = 5 in a finite box of size L = 8. The simulation was performed using No = 4, 
with couplings Cq = 0.487259, C2 = 0.298043, C4 = -0.211675 and Cg = 0.0405311. The 
temporal extent of the lattice was chosen to be T = 64, which is approximately three times 
larger than the predicted inverse energy difference in the ground and first excited states 
and A2, given by Liischer's formula. The lowest two measured energies, Ai and A2, were then 



determined by fitting the time-dependence of the effective mass (Eq. CI with At = 1) to 
a constant plus exponential form. The effective mass was fit over an interval [Tmin,'''max], 
where Tmax = T and r^m was varied until a plateau was achieved in the fit values of Ai and 
A2. We found that the results for the ground state and first excited state agreed with the 
theoretical result determined by Liischer's formula to within errors of about 0.007% and 8% 
percent, respectively. 

The implementation of the external potential was checked by measuring the ground and 
excited state energies of a single fermion confined to a harmonic trap as a function of L/ Lq. 
In this case, since there are no auxiliary fields present, this check may be regarded as a 
numerical calculation, rather than a simulation. We chose Lq = 2,2.5,3,3.5 and 4, and 
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oj = 0.005 in order to ensure negligible temporal and spatial discretization errors. Sources 
with appropriate symmetry properties were constructed in order to extract the lowest seven 
energies of the SHO. In all calculations, the temporal extant was chosen such that T ^ l/u in 
order to eliminate contamination from higher energy states in the single fermion correlation 
functions. Numerical calculations obtained from spatial lattices L = 8, 16, 32, 64 are shown 
in Fig. [5} and show good agreement with the continuum theory. 



[1] H. Hu, P. D. Drummond, and X.-J. Liu, Nature Physics 3, 469 (2007), arXivxond- 
mat/0701744. 

[2] D. B. Kaplan, M. J. Savage, and M. B. Wise, Nucl.Phys. B534, 329 (1998), nucl-th/9802075. 
[3] D. B. Kaplan, M. J. Savage, and M. B. Wise, Phys.Lett. B424, 390 (1998), nucl-th/9801034. 
[4] K. M. O'Hara, S. L. Hemmer, M. E. Gehm, S. R. Granade, and J. E. Thomas, Science 



298, 2179 (2002), http://www.sciencemag.org/content/298/5601/2179.full.pdf, URL [http: 
|//www . sciencemag . org/ content/298/5601/2179 . abstract! 
[5] C. A. Regal and D. S. Jin, Phys. Rev. Lett. 90, 230404 (2003). 

[6] M. E. Gehm, S. L. Hemmer, S. R. Granade, K. M. O'Hara, and J. E. Thomas, Phys. Rev. A 

68, 011401 (2003). 
[7] T. Bourdel et al., Phys. Rev. Lett. 91, 020402 (2003). 

[8] K. Dieckmann, C. A. Stan, S. Gupta, Z. Hadzibabic, C. H. Schunck, and W. Ketterle, Phys. 

Rev. Lett. 89, 203201 (2002). 
[9] C. A. Regal, M. Greiner, and D. S. Jin, Phys. Rev. Lett. 92, 040403 (2004). 
[10] J. Kinast, S. L. Hemmer, M. E. Gehm, A. Turlapov, and J. E. Thomas, Phys. Rev. Lett. 92, 

150402 (2004). 

[11] C. A. Regal, C. Ticknor, J. L. Bohn, and D. S. Jin, Nature (London) 424, 47 (2003), 

arXiv:cond-mat/0305028. 
[12] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. Raupach, A. J. Kerman, and W. Ketterle, 

Phys. Rev. Lett. 92, 120403 (2004). 
[13] M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. H. Denschlag, and R. Grimm, 

Phys. Rev. Lett. 92, 120401 (2004). 
[14] P. Arnold, J. E. Drut, and D. T. Son, Phys. Rev. A 75, 043605 (2007), arXivxond- 



46 



mat/0608477. 

[15] S. Tan, ArXiv Condensed Matter e-prints (2004), arXiv:cond-mat/0412764. 

[16] Y. Nishida and D. T. Son, Physical Review Letters 97, 050403 (2006), arXiv:cond- 

mat/0604500. 

[17] Y. Nishida and D. T. Son, Phys. Rev. D 76, 086004 (2007), 0706.3746. 
[18] F. Werner and Y. Castin, Phys. Rev. A 74, 053604 (2006), arXiv:cond-mat/0607821. 
[19] D. T. Son and M. Wingate, Annals of Physics 321, 197 (2006), arXiv:cond-mat/0509786. 
[20] E. Braaten, D. Kang, and L. Platter, Phys.Rev. A78, 053606 (2008), 0806.2277. 
[21] D. Lcc and T. Schafcr, Phys. Rev. C 72, 024006 (2005), arXiv:nucl-th/0412002. 
[22] A. Bulgac, J. E. Drut, and P. Magierski, Phys.Rev.Lett. 96, 090404 (2006), cond-mat/0505374. 
[23] D. Lee and T. Schafer, Phys.Rev. C73, 015201 (2006), nucl-th/0509017. 
[24] D. Lee and T. Schafer, Phys.Rev. C73, 015202 (2006), nucl-th/0509018. 
[25] D. Lee, Phys. Rev. C78, 024001 (2008), 0803.1280. 
[26] A. Gezerlis and J. Carlson, Phys. Rev. C 77, 032801 (2008), 0711.3006. 
[27] A. Bulgac, J. E. Drut, and P. Magierski, Phys. Rev. A 78, 023625 (2008), 0803.3238. 
[28] M. McNeil Forbes, S. Gandolfi, and A. Gezerlis, ArXiv e-prints (2010), 1011.2197. 
[29] L. Pricoupenko and Y. Castin, Journal of Physics A Mathematical General 40, 12863 (2007), 
0705.1502. 

[30] G. E. Astrakharchik, J. Boronat, J. Casulleras, and S. Giorgini, Physical Review Letters 93, 

200404 (2004), arXiv:cond-mat/0406113. 
[31] S. Gandolfi, K. E. Schmidt, and J. Carlson, Phys. Rev. A 83, 041601 (2011), 1012.4417. 
[32] J. Carlson, S.-Y. Chang, V. R. Pandharipande, and K. E. Schmidt, Phys. Rev. Lett. 91, 

050401 (2003). 

[33] O. Goulko and M. Wingate, Phys. Rev. A 82, 053621 (2010), 1008.3348. 

[34] E. Burovski, N. Prokof'ev, B. Svistunov, and M. Troyer, Phys. Rev. Lett. 96, 160402 (2006). 

[35] A. Bulgac, Phys. Rev. A 76, 040502 (2007), arXiv:cond-mat/0703526. 

[36] T. Papenbrock, Phys. Rev. A 72, 041603 (2005), arXiv:cond-mat/0507183. 

[37] S. Bour, X. Li, D. Lee, U. Meifiner, and L. Mitas, ArXiv e-prints (2011), 1104.2102. 

[38] T. Abe and R. Seki, Phys.Rev. C79, 054003 (2009), 0708.2524. 

[39] J. Carlson and S. Reddy, Phys. Rev. Lett. 95, 060401 (2005). 

[40] A. Gezerlis and J. Carlson, Phys.Rev. C77, 032801 (2008), 0711.3006. 

47 



[41] D. Blume, J. von Stecher, and C. H. Greene, Phys. Rev. Lett. 99, 233201 (2007). 
[42] S. Chang and G. Bertsch, Phys.Rev. A76, 021603 (2007), physics/0703190. 
[43] S. Chang, V. Pandharipande, J. Carlson, and K. Schmidt, Phys.Rev. A70, 043602 (2004), 
physics/0404115. 

[44] M. G. Endres, D. B. Kaplan, J.-W. Lee, and A. N. Nicholson, PoS LATTICE2010, 182 
(2010), 1011.3089. 

[45] J.-W. Lee, M. G. Endres, D. B. Kaplan, and A. N. Nicholson, PoS LATTICE2010, 197 
(2010), 1011.3026. 

[46] A. N. Nicholson, M. G. Endres, D. B. Kaplan, and J.-W. Lee, PoS LATTICE2010, 206 
(2010), 1011.2804. 

[47] J.-W. Chen and D. B. Kaplan, Phys.Rev.Lett. 92, 257002 (2004), hep-lat/0308016. 
[48] M. G. Endres, D. B. Kaplan, J.-W. Lee, and A. N. Nicholson (2011), 1106.0073. 
[49] M. G. Endres, D. B. Kaplan, J.-W. Lee, and A. N. Nicholson (to appear). 
[50] M. Liischer, Commun.Math.Phys. 104, 177 (1986). 
[51] M. Liischer, Commun.Math.Phys. 105, 153 (1986). 
[52] M. Liischer, Nucl.Phys. B354, 531 (1991). 

[53] S. Beane, P. Bedaque, A. Parreno, and M. Savage, Phys.Lett. B585, 106 (2004), hep- 
lat/0312004. 

[54] M. C. Birse, Phil.Trans.Roy.Soc.Lond. (2010), 1012.4914. 

[55] H. W. Griesshammer, Nucl.Phys. A760, 110 (2005), nucl-th/0502039. 

[56] H. W. Griesshammer, Few Body Syst. 38, 67 (2006), nucl-th/0511039. 

[57] D. Blume and K. Daily, Comptes Rendus Physique 12, 86 (2011), ISSN 1631-0705, few body 
problem, URL |http : //www . sciencedirect . com/science/article/B6X19-51XNXCG- 1/2/ 
Ibl72a44ble4a2b250bae52fcfeca2f31[ 

[58] D. Lee, Eur.Phys.J. A35, 171 (2008), 0704.3439. 

[59] M. Liischer, Comput.Phys.Commun. 79, 100 (1994), hep-lat/9309020. 



48 



